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Abstract 

We use molecular dynamics (MD) to simulate an unstable homogeneous mixture of binary fluids 
(AB), confined in a slit pore of width D. The pore walls are assumed to be flat and structureless, 
and attract one component of the mixture (A) with the same strength. The pair-wise interactions 
between the particles is modeled by the Lennard-Jones potential, with symmetric parameters that 
lead to a miscibility gap in the bulk. In the thin-film geometry, an interesting interplay occurs 
between surface enrichment and phase separation. 

We study the evolution of a mixture with equal amounts of A and B, which is rendered unstable 
by a temperature quench. We find that A-rich surface enrichment layers form quickly during the 
early stages of the evolution, causing a depletion of A in the inner regions of the film. These 
surface-directed concentration profiles propagate from the walls towards the center of the film, 
resulting in a transient layered structure. This layered state breaks up into a columnar state, 
which is characterized by the lateral coarsening of cylindrical domains. The qualitative features 
of this process resemble results from previous studies of diffusive Ginzburg-Landau-type models 
[S. K. Das, S. Puri, J. Horbach, and K. Binder, Phys. Rev. E 72, 061603 (2005)], but quantitative 
aspects differ markedly. The relation to spinodal decomposition in a strictly 2-d geometry is also 
discussed. 

PACS numbers: 68.05.-n,64.75.+g,68.08.Bc,68.15.+e 
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I. INTRODUCTION 



Thin fluid films have a broad range of applications in technology as lubricants, protecting 
layers, for production processes of layered structures in microelectronics, etc. In particular, 
ultra-thin films have become extremely important in the context of nanotechnolog y 1 ^ 3 ! 4 ^ 6 ! 7 . 
The interplay of surface effects and finite-size effects with the bulk behavior of these systems 
poses challenging theoretical problem o 8 ! 9 ' 10 ' 11 . A particularly interesting problem in this 
context is the phase separation of binary (or multi-component) mixtures in thin films, as 
most materials of practical interest have more than one component, e.g., metallic alloys, 
ceramics, polymer blends, etc. Phase changes in reduced geometries (e.g., 2-d systems) differ 
in many aspects from the bulk behavior in three dimensions. The interplay between surface 
and bulk behavior leads to complex phenomena such as wetting transitions, prewetting 
and layering transitions, etc . 12 i 13 i 14 i 15 i 16 i 17 i 18 i 19 i 20 . These transitions may compete with phase 
changes that occur in the bulk, such as phase separation in mixture a 21 i 22 i 23 i 24 i 25 i 26 . There has 
been intense study of phenomena such as surface- directed spinodal decomposition (SDSD) 
or surface- directed phase 5epara^ori i 27 i 28 i 29 i 30 i 31 i 32 i 33 i 34 i 35 , but our theoretical understanding of 
these problems is still incomplete. 

To gain a better understanding of SDSD in thin fluid films, we have undertaken a com- 
prehensive molecular dynamics (MD) simulation of binary mixtures in a slit pore. A pre- 
liminary account of our results has been published as a letter—. In this paper, we present 
detailed results from this study. Our MD simulations are based on a symmetric binary (AB) 
Lennard- Jones (LJ) mixture, for which the bulk phase diagram has been determined to high 
accuracy 3 ^. We also have a good understanding of various other properties of this mix- 
ture, e.g., static and dynamic response and correlation functions, transport coefficient o 37 ' 38 , 
and the interfacial tension between coexisting A- rich and B-rich phases 39 . Many previ- 
ous simulation studies of phase-separation kinetics in a thin-film geometr y 40 ! 41 have used 
Ginzburg-Landau (GL) models, and hence lack any direct connection to a microscopic de- 
scription. Our present modeling bridges the gap between the atomistic description of liquids 
and the mesoscale domain structures that form as the kinetics of phase separation proceeds 
(see Fig. Q). This direct approach has the further merit that hydrodynamic interactions 
are automatically incorporated. It is well-known that these interactions have a pronounced 
effect on the kinetics of domain growt h 21 ^ 22 ! 23 ! 24 ! 25 ! 26 ! 42 ! 43 ! 44 ! 45 ! 46 . 
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In order to account for hydrodynamics in the framework of GL models, rather extensive 
computations are require d 47 ! 48 ' 49 ' 50 ' 51 . The GL approach is appropriate if one is primarily 
concerned with the scaling behavior of the late-stage domain growt h 21 ' 22 ' 23 ' 24 ' 25 ' 26 . But re- 
cent work based on lattice Boltzmann simulations questions the quantitative validity of GL 
models, and reports very slow crossovers extending over many decades in time 5 ^ (however, 
a "somewhat narrower" crossover region was reported in a subsequent work^ 3 -). This study 
also raises questions about some of the previous work on this problem, using the lattice 
Boltzmann method or related "lattice gas" -approximations to the Navier-Stokes equations 
of hydrodynamic o 54 ' 55 ' 56 ' 57 ' 58 . However, the lattice Boltzmann description is even more re- 
mote from an atomistic description of matter than the GL model. Further, it is not clear 
how one incorporates the proper boundary conditions with respect to complete or partial 
wetting in such an approac h 59 ' 60 . The studies mentioned above are primarily concerned with 
domain growth in d = 3. For the d = 2 case, it is even controversia l 61 ' 62 ' 63 to what extent 
a scaling behavior describes the late stages of coarsening. In view of these problems, it is 
useful to undertake an MD simulation despite the fact that the accessible scales in length 
and time are limited. Earlier MD studies of bulk phase separatio n 64 ' 65 have addressed coars- 
ening in d = 3 and in d — 2, but the latter work is somewhat inconclusive^. Further, there 
have only been preliminary MD studies of SDSD in thin film o 67 ' 68 . 

The outline of this paper is as follows. In Sec. II, the theoretical background (equilibrium 
phase behavior of binary mixtures in a thin-film geometry, theory of domain growth, etc.) 
will be concisely reviewed. In Sec. Ill, we provide details of our MD methods. In Sec. IV, 
we present simulation results for SDSD in thin films. Finally, Sec. V concludes this paper 
with a summary and discussion, including a comparison to the GL approach. 

II. THEORETICAL BACKGROUND 

A. Equilibrium phase behavior of binary mixtures confined between walls 

A homogeneous binary mixture becomes unstable to phase separation when it is quenched 
into the miscibility gap (see Fig. |2j). For a symmetric mixture, the miscibility gap is sym- 
metric with respect to the concentration x™ 1 * = 1/2. 

At the surface of a semi-infinite mixture, one may encounter a wetting 
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transitio n 12 i 13 i 14 i 15 i 16 i 17 i 18 i 19 i 20 . This transition implies a singular behavior of the surface 
excess free energy F$, which is defined as (for a film between two walls at distance D) 
Fftim = F h + 2F S /D, D — > oo, F h being the bulk free energy of the system. Assuming, as 
done in Fig. |21 that the wetting transition occurs at the surface of B-rich mixtures (caused 
by the preferential attraction of A-particles to the walls), the transition is characterized by a 
divergence of the surface excess concentration of A, X s ™ . This quantity can be obtained from 
F$ via suitable derivatives, or by integrating the concentration profil o 12 i 13 i 14 i 15 i 16 i 17 i 18 i 19 i 20 i 69 

D/2 

< rf = / [xAtf-xW^dz , D^OC , (1) 



where z is the distance from the wall, which is located at z — 0. If the wall is nonwet (or 
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partially we£) 12 i 13 i 14 i 15 i 16 i 17 i 18 i 19 i 20 , x™ rf tends to a finite value (^a cLx) when xa — > x 
from the one-phase region. On the other hand, for a wet (or completely wet) wall, x™ rf = oo 
- corresponding to an infinitely thick A-rich wetting layer coating the wall, separated from 
the B-rich bulk by a flat interface. 

At the coexistence curve x^ coex , the surface excess free energy is that of an A-rich phase 
^Scoex* if the W& 11 is nonwet. For a wet wall, we have F$ = F^~J^ + 7ab, 7ab being 
the interfacial tension between coexisting A-rich and B-rich phases. These quantities also 
determine the contact angle 9 at which an A-B interface in the nonwet region meets the 

-■-J2l2.13.14.1B.1B.17.1H.19.2n 

/ pB- rich z^A— rich\ 

„„„ a _ V S,coex S.coex ) . f pB-rich , pA-rich , /n\ 

COSU- , It ^ s < ^ S ,coex +7AB • {*) 

7AB 

If the state of the system is changed such that one increases the temperature but stays 
always at the coexistence curve a^loex' one encounters a wetting transition at temperature 
T w (Fig. I2J), where the state of the wall changes from nonwet (T < T w ) to wet (T > T w ). This 
transition may be of second order (Fig. |2t>) or first order (Fig. Eb) • in the second-order case, 
x surf ,-[i ver g es continuously when T — > T~, while otherwise there is a discontinuous jump in 
^surf f rom a f| m te value at T w to oo at T+. In the first-order case, there is also a prewetting 
transition in the one-phase region (Fig. db), where the thickness of the A-rich surface layer 
jumps from a smaller value to a larger (but finite) value. This line of prewetting transitions 
ends in a prewetting critical point. 

This brief review of wetting phenomena provides the basis to understand the equilibrium 
behavior of binary mixtures in thin film o 70 ! 71 ! 72 ! 73 ! 74 ! 75 ! 76 ! 77 ! 78 ! 79 . If the walls are neutral (i.e., 



it has the same attractive interactions with both A-particles and B-particles), the critical 
concentration remains x c A rit = 1/2. However, the critical temperature T C (D) is lowere d 70 ' 75 ' 76 
relative to the bulk: 

T c - T C (D) oc D~ x '\ (3) 

where v ~ 0.62 9 80 i 81 is the critical exponent of the correlation length £ of concentration 
fluctuations (in the universality class of the d = 3 Ising model). Note, however, that critical 
correlations at fixed finite D can become arbitrarily long-range only in the lateral direction 
parallel to the film. Thus, the transition at T C (D) belongs to the class of the d = 2 Ising 
model. The states below the coexistence curve of the thin film correspond to two-phase 
equilibria characterized by lateral phase separation. 

When there is a preferential attraction of A-particles to the walls, the phase diagram of 
the thin film is no longer symmetric with respect to xa = 1/2, although we did assume such 
a symmetry in the bulk. The shift of x c A Tlt and the resulting change of the coexistence curve, 
is the analog of capillary condensation of gase o 73 ^ 2 for binary mixtures. 

The coexisting phases in the region below the coexistence curve of the thin film are 
inhomogeneous in the direction perpendicular to the walls (see Fig. |2t). In the A-rich 
phase, we expect only a slight enhancement of the order parameter ip(z), which is defined 
in terms of the densities ua(z), Ub(z) of A and B particles as 

n A (.)-n B (z) 
n A (z) +n B (z) 

In the B-rich phase, however, we expect pronounced enrichment layers. As D — > oo, the 
thickness of these layers diverges for T > T w but stays finite for T < T w . In a film of finite 
thickness, the width of A-rich surface layers also stays finite, e.g., for T > T w , x A oc InD 
for short-range surface forces, while x s A ri oc D 1 ^ 3 for non-retarded van der Waals' force d 75 ' 8 ? . 
Thus, the wetting transition is always rounded off in a thin film. The prewetting line (Fig. 12b) 
does have an analog in films of finite thickness D, for sufficiently large D. This transition 
splits into a two-phase region at small xa between the thin-film triple point and the thin-film 
critical point on the B-rich side. This two-phase region corresponds to a coexistence between 
B-rich phases with A-rich surface layers, both of which have finite (but different) thickness. 
As D — > oo, the thin-film critical point on the B-rich side moves into the prewetting critical 
point, while the thin-film triple point merges with the first-order wetting transition. On the 
other hand, when D becomes small, the thin-film critical point and the thin-film triple point 



may merge and annihilate each other. For still smaller D, the thin-film phase diagram then 
has the shape shown in Fig. |2t, although one has first-order wetting in the semi-infinite bulk 
(Fig. 13)). 

Finally, we comment on the state encountered below the bulk coexistence curve, but 
above the coexistence curve of the thin film. When one crosses the bulk coexistence curve, 
there is a rounded transition towards a layered (stratified) structure with two A-rich layers at 
the walls and a B-rich layer in the middle. The temperature range over which this rounded 
transition is smeared is also of order AT oc D~ x l v around T c . Hence, for large D, this 
segregation in the direction normal to the walls may easily be mistaken (in experiments or 
simulations) as a true (sharp) phase transition. We stress that this is not a true transition 
- one is still in the one-phase region of the thin film, although the structure is strongly 
inhomogeneous! The situation qualitatively looks like the concentration profile shown in the 
upper part of Fig. |2t- The difference is that, for D — > oo, the thickness of true wetting 
layers scales sub-linearly with D, as noted above. However, for phase separation in the 
normal direction which gradually sets in when one crosses the bulk coexistence curve, one 
simply has A-rich domains of macroscopic dimensions (proportional to D) adjacent to both 
walls. Unfortunately, the layers resulting in this stratified structure are often referred to as 
"wetting layers" in the literature, although this is completely misleading. We reiterate that 
A-rich wetting layers only form when a B-rich domain extends to the surface, which is not 
the case here. 

We also caution the reader that a picture in terms of A-rich layers at the walls and a B-rich 
domain in the inside of the film is an over-simplification because the thickness of the domain 
walls cannot really be neglected in the region T C (D) < T < T c , where a stratified structure 
occurs in equilibrium. This is seen from the relation £ oc (1 — T/T c )~ u , in conjunction with 
Eq. Q, which shows that £ ~ 0(D) at T C (D). Thus, domains and domain walls are not 
well-distinguished in the region under consideration, since the interfacial width is 0(£) 20 ' 69 . 

When the interface between A-rich and B-rich domains is treated as a sharp kink (this 
approximation is popular in theoretical treatments of wettin g 12 ! 1 ^ 14 ' 1 ^ 1 ^ 1 ^ 1 ^ 19 ' 21 ^ ), one might 
think that a sharp wetting transition could still be described in terms of the vanishing of 
the contact angle 9 as T — > T w (Fig. Efc). However, it is clear that for a correct treatment 
the finite width of the interface needs to be taken into account. Thus, for finite D, the 
contact angle in Fig. Efc is ill-defined, and the transition between the two states depicted 
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in Fig. is smooth, because a B-rich nonwet domain may also have a thin A-rich layer at 
its surface (x™ rf , in general, is nonzero). One should also note that the contact "line" is 
distorted by line tension effects when it hits the wall, and the line tension of the interface 
at the wall would also modify Eq. (J2)| 84 i 85 i 86 . The difficulty of estimating the contact angle 
in finite geometries is well-known from studies of nanoscopic droplet o 87 i 88 . 

The central conclusion in this subsection is that in the final equilibrium to which, for 
times t — > oo and for small D, the thin film evolves, there is no fundamental difference 
whether or not we are above or below the wetting transition temperature, but it matters 
whether T < T C (D) or T > T C {D). 



B. Bulk phase separation of binary fluid mixtures 

Next, we review our understanding of the kinetics of phase separation in bulk fluid mix- 
tures which are rendered thermodynamically unstable by a rapid quench (at t = 0) into the 
miscibility gap (see Fig. EJ). The initial state (t < 0) is spatially homogeneous, apart from 
small-scale concentration inhomogeneities. The final equilibrium state consists of macro- 
scopic domains of the two coexisting phases, with relative amounts determined by the 
lever rule. We are interested in the evolution from the initial homogeneous state to the 
final segregated state. For quenches below the spinodal curve, the homogeneous system 
is unstable and decomposes via the spontaneous growth of long-wavelength concentration 
fluctuation o 21 i 22 i 23 i 24 i 25 i 26 (spinodal decomposition). Understanding the full time evolution 
from the initial stages to the late stages of coarsening is a formidable problem, and is typi- 
cally accessed by large-scale simulations of coarse-grained models. 

Nevertheless, there exist some cases in which simple domain growth laws can be obtained 
from analytical consideration ;] 4 ? i44i4 ^ 4 ^ 89i9 ° . The evaporation- condensation mechanism of 
Lifshitz and Slyozov (LS) 89 corresponds to a situation where a population of droplets of the 
minority phase (say, A) is in local equilibrium with the surrounding supersaturated majority 
phase. The LS mechanism leads to a growth law (valid for dimensionality d > 1) £{t) oc t 1 ^ 3 , 
t — > oo, where £(t) is the linear dimension of the droplets. 

The droplet diffusion- coagulation mechanism^ is specific to fluid mixtures, and is based 
on Stokes law for the diffusion of droplets, yielding 9 - £(t) oc (t/r)) 1 ' d . 

A faster mechanism of domain growth in fluids was proposed by Siggia 4 ^, who studied 
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the coarsening of interconnected domain structures via the deformation and break-up of 
tube-like regions, considering a balance between the surface energy density ~ 7ab/-^ and the 
viscous stress ~ 67ir]Ve/£ 2 ^. Thus, vi oc Jab/v and ^ cx vg, or £{t) oc 1 ^-t in d = 3. In d = 2, 
the analog of this hydrodynamic mechanism is controversial. San Miguel et al^ argue that 
strips (d = 2 analogs of tubes) are stable under small perturbations, in contrast to the d = 3 
case. For critical volume fractions, an interface diffusion mechanism is proposed which yields 
£{t) oc t 1/2 , i.e., the same growth law as the Brownian coalescence mechanism of droplets in 
d = 2 (see above). On the other hand, Furukaw a 45 i 62 argues for a linear relation £(t) oc t in 
d = 2 as well. However, recently there is growing evidenc e^ 1 that different characteristic 
length scales in d = 2 may exhibit different growth exponents, suggesting that there is no 
simple dynamical scaling of domain growth in d = 2! 

Finally, we remark that the above growth laws do not constitute the true asymptotic 
behavior, either in d = 2 or d = 3. Rather, these results only hold for low enough Reynolds 
numbers^. For £ £i n — ?7 2 /( n 7AB), the so-called inertial length 2 ^, one enters a regime 
where the surface energy density 7ab/^ is balanced by the kinetic energy density nvf. This 
yields the following growth law for the inertial regim e 2 ^ : 



which is valid for both d = 2 and d = 3. In d = 2, evidence for both £{t) oc t 1 ^ 2 and 
£(t) oc t 2//3 has been reported, but the conditions under which such power laws hold in d = 2 



C. Ginzburg-Landau model of surface-directed spinodal decomposition 

In this subsection, we briefly discuss a coarse-grained description of binary mixtures 
in a thin-film geometry, which can reproduce the phase diagrams shown in Fig. This 
description can also be used to obtain a model for the kinetics of phase separation in a 
confined geometry. Our starting point is a mean-field description of a binary mixture near 
its critical point, where one introduces a local order parameter ip(p,z). (Here, p represents 
the coordinates parallel to the walls, and z is the coordinate in the perpendicular direction, 
as before.) The surfaces of the thin film S\ and S2 are located at z = and z = D, 
respectively. 




(5) 



are still not clear 51 i 61 i 62 i 63 ' 64 i 65 ' 66 . 
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We denote the order parameter describing the bulk coexistence curve as ifo, and define 
i/j'(p,z) = ip(p, z)/ipb. Further, we measure distances p,z,D in units of £. Then, the 
dimensionless free-energy functional of a binary mixture in a thin-film geometry can be 
written as a sum of a bulk term and two surface terms^i, F [ip] = F^ [ip] + Fg 1 [ip] + F$ 2 [ijj] , 
where we have dropped the prime on if>'. Here, we have 

D 

dp / dz 



o 



2 4 4 



(6) 



The terms Fg 1 and F$ 2 are obtained as integrals over the surfaces Si and S2: 
F 3l = f rfp{-|[^(p,0)] 2 -^(p,0)- 7 ^(p,0)^ 



Si 



z=0 



(7) 



and analogously for Fs 2 . 

In Eq. (jHJ), we have included a ^-dependent surface potential which arises due to the 
surfaces. In our subsequent discussion, we will consider symmetric power-law potentials: 

V(z) = -V [(z + 1)-p+(D + 1-z)~p] , (8) 

which satisfy V(z) = V(D — z). The potentials are taken to originate behind the surfaces 
so as to avoid singularities at z = 0, D. 

The terms Fs iy Fs 2 represent the surface excess free-energy contributions due to local 
effects at the walls, with g and 7 phenomenological paramete ra 31 ' 35 ' 69 ' 91 . The dimensionless 
surface fields in Fs 1 and F$ 2 are hi = — V(0) and h 2 = —V(D), respectively. The one-sided 
derivatives appear in Fs 1 and F$ 2 due to the absence of neighboring atoms for z < and 
z > D. 

Let us first consider the limit D — > 00. For the long-range surface potential Eq. (jSJ), 
only first-order wetting transitions are possible^. For power-law potentials as in Eq. (jHJ), an 
approximate theory 3 ^ predicts that the wetting transition occurs when 2Vq/ (p — 1) = 7ab- 

For finite values of D, one can obtain phase diagrams of thin films (as shown schemat- 
ically in Fig. 12} by minimizing the free-energy functional in Eqs. (E}-(JH}- However, this 
requires numerical wor k 70 1 77 . We also note that the ?/> 4 -Hiodel cannot describe either the 
low-temperature region (where complete separation between A and B occurs), or the non- 
mean-field critical behavior. 
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We next discuss the dynamics of phase separation in thin films. First, let us estab- 
lish the dynamical equations which govern phase separation in the bulk for the diffusion- 
driven case. The local order parameter ip(r,t) is conserved, and obeys the continuity 



equation 21 i 22 i 23 i 24 i 25 i 26 : 



^il>(f,t) = -V • J{r,t) , (9) 



The current J(r, t) contains contributions from the local chemical potential difference /i(f, t), 
and from statistical fluctuations, 0(f,t): 

J(r,t) = -Vfi{r,t) + 6{r,t) , 

= sWT) • (10) 

Using the ^ 4 -free-energy functional in Eq. (JB|) . we obtain the dynamical model: 

9(f,t)\, 0<z<D. (11) 



^(f,t) = V-<j V 



-ip + ip 3 - -V 2 ^ + V(z) 
2 



We assume that the noise 9 is a Gaussian white noise, obeying the relations (9(r,t)) = 
and (9i(f ',t')9j(f ",<")) = 2e5 ij 5{f ' - f ")5(f - t"), where the indices i,j denote the 
Cartesian components of vector 9. Note that the time units have been chosen such that the 
diffusion constant in Eq. (jlUj) is unity. With respect to the dynamical behavior in the critical 
region, Eq. (|11|) corresponds to model B in the Hohenberg-Halperin classifications 2 .. However, 
statistical fluctuations are irrelevant for the late stages of spinodal decomposition 93 : The 
deterministic model obtained by setting e = in Eqs. (J5 )) -([1U) ) also yields the LS growth law 
(see Sec. IIB) in the late stages of domain growth. 

One can incorporate hydrodynamic effects, as is appropriate for fluid mixtures, by in- 
cluding a velocity field v(r, t)^, but this will not be further considered here. 

The above models describe coarsening kinetics in the bulk. When one deals with SDSD in 
thin films, the model needs to be supplemented by boundary conditions at the surface d 28 ! 94 . 
The first boundary condition expresses the physical requirement that the z-component of 
the flux at the surfaces must vanish: 



J z (p,0,t) 



d_ 

dz 



-tfj + tfj 3 - -Vfy + V(z) 



+ 9A =0 , (12) 



2=0 



and similarly for z = D. The second boundary condition describes the evolution of the 
surface order parameter. Since this quantity is not conserved, it is described by a relaxational 
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kinetics of model A type**: 



9 ,,^ n s SF 



<9t r vr ' ' ' 5ijj(p,0,t) 

= hi+giJ>(p,0,t)+7^ ■ (13) 

OZ 2=0 

An analogous equation can be written down for the relaxation of ip(p, D, t). Here, To sets the 
time-scale of this nonconserved kinetics. Since ip{p, 0, t) relaxes much faster than the order 
parameter in the bulk, it is reasonable to set tq = C 35 . Then, the dynamics as well as the stat- 
ics is controlled by two surface parameters, and g/j. During the early stages of SDSD, 
the fast relaxation of the order parameter at the surfaces provides a boundary condition for 
the phase of concentration waves that grow in the thin film. In the bulk, the random orien- 
tations and phases of these growing waves do not yield a systematic evolution of the average 
order parameter. However, the surface-directed concentration waves add up to give an aver- 
age oscillatory concentration profile near the surfaces of a thin fil m 27 i 28 i 2 V°i 31 i 32 i 33 i 34 i 35 i 40 i 41 . 

Our early work on SDSD in thin films 40 omitted both hydrodynamic interactions and the 
noise in Eq. (J1U)) . and focused on the d = 2 case. In recent work 41 -, we have studied the 
d = 3 case using the GL model in Eq. (fTTj) with the noise term, in conjunction with the 
boundary conditions in Eqs. (fT2"j) - (fT3~|) . In Sec. IV, we will compare our MD results with 
results from this study. The details of the GL simulation are as follows. We implemented an 
Euler-discretized version of Eqs. (111}) . (jl2j) - (fT3*|) on an L x L x D lattice. The discretization 
mesh sizes were Ax = 1 and At = 0.02. The surface potential was of the form in Eq. (JHJ) with 
p = 3, which corresponds to a non-retarded van der Waals' interaction between the surfaces 
and a particle in d = 3. The parameter values were g = —0.4,7 = 0.4, and Vq = 0.325 for 
D = 5 and Vq = 0.11 for D = 10, corresponding to a partially wet surface in equilibrium. 
We stress that, for a fluid mixture, the above diffusive model is relevant during the early 
stages of phase separatio n 21 i 22 i 23 i 24 i 25 i 26 , but is not expected to yield useful results for the 
intermediate and late stages of domain growth. 



III. MODEL AND MOLECULAR DYNAMICS METHODS 

For our MD study, we consider a fluid of point particles located in continuous space in a 
box of volume L x L x D. We apply periodic boundary conditions in the x and y directions, 
while impenetrable walls are present at z = and z = D. These walls give rise to an 
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integrated LJ potential (a = A,B): 



u w (z) 



2%na 3 
3 



2 /a\9 /a\3' 
— o, 



(14) 



15 Vz'/ " \z'. 

where n is the reference density of the corresponding bulk flui d 37 i 38 i 39 , a is the LJ diameter 
of the particles, and e w is an energy scale for the strength of the wall potentials. Further, 
5a = 1 and 5b = 0, so A-particles are attracted by the walls while B-particles are not. The 
coordinate z' = z + a/2 for the wall at z = 0, and z' = D + cr/2 — z for the wall at z = D. 
Therefore, the singularities of u w (z) do not occur within the range < z < D but rather at 
z = —a/2 and z = D + a/2, respectively. 

The particles in the system interact with LJ potentials: 



ufa) = 4e Q/3 



12 / \ 6 
Gap \ / 



ra = \ri-rj\, (15) 



where a, (3 = A, B. The LJ-parameters (e Qj a, a a p) are chosen as follows: 

oaa = aAB = 0"BB = 0", 

£aa = cbb = e, e AB = ~ • (16) 

The units of length, temperature, and energy are chosen such that a — 1, e = 1, hs — 1- 
The masses of the particles are chosen to be equal, m A = m-Q = m = 1. Thus, the MD time 
uni1i"2i: 

fw 2 \ 1/2 1 

to H^J (17) 

becomes a dimensionless number. To speed up the calculations, the LJ potential is truncated 
at Tij = 2.5 a and shifted to zero there, as usual^. To ensure that our study of fluid-fluid 
phase separation is not affected by other phase transitions (e.g., liquid-gas or liquid-solid 
transitions), we work with bulk density n = 1.0 and focus on temperatures T > 1.0. In 
principle, in thin films one could have a wall-induced crystallization at temperatures above 
the bulk melting temperature: however, we have not seen any evidence for such an effect in 
our model. In our previous work on the bulk behavior of the same mode l 3 ^ 3 ^ 3 ? , we found 
that the critical temperature for bulk phase separation is T c ~ 1.638. Here, we present results 
from simulations of quenching experiments to T = 1.1. At this temperature, the bulk phase 
separation is essentially complete. In addition, the bulk correlation length £ ~ 1 (i.e., one LJ 
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diameter) within the relative accuracy of about 5% to which it can be determined^. Further, 
material parameters which enter the theories reviewed in Sec. II (e.g., the interfacial tension 
7ab, the shear viscosity rf) are explicitly known as well. The appropriate values are jab — 0.9 
(see Ref.—) and rj ~ 7 (see Ref.—). (Recall that these quantities are measured in LJ units 
and hence are dimensionless.) The availability of most material parameters for our system is 
a distinct advantage of our atomistic model in comparison to coarse-grained models, where 
it is often unclear what ranges of effective parameters correspond to physically reasonable 
choices. 

The strength of the wall-particle interaction is taken as e w = 0.005, which corresponds to 
partially wet walls at T = 1.1. To find the precise location of the wetting transition for our 
model would require a major computational effort, and this has not been attempted. Recall 
that one expects T w — ► T c when e w — > - this was the motivation for choosing a rather 
small value of e w in our study. However, due to the special choices made [such as Eq. ffTBj) ]. 
for the sake of simplicity, it would be premature to try to explicitly relate our model to a 
specific real system. 

The lateral size L of the simulated systems must be large enough that the laterally 
inhomogeneous structures that form during segregation are not affected by finite-size effects. 
Therefore, we chose L = 128 for the thinnest film in our study (D = 5) and L = 64 for the 
thicker ones (D = 10, 20). As the confining potentials diverge at z — — cr/2 and z = D-\-a/2 
(a = 1), the volume in which particles can be is V = L 2 (D + 1). We will report results 
from three sets of simulations, with D = 5, L = 128 (N = 98304 particles); D = 10, L = 64 
(N = 45056); and D = 20, L = 64 (N = 86016). Thus, the particle density is n = 1 in 
all these cases. For the range of times studied here (t < 8000), test runs with other linear 
dimensions showed that our choices of L are large enough to eliminate finite size effects, 
within the limits of our statistical accuracy. Of course, for a study on larger time scales also 
larger system sizes would be required! 

The initial states of the simulations need to be carefully prepared. We equilibrated a fluid 
of iV particles (with TVa = ^Vb = N/2) in the specified volume at a very high temperature 
(T = 5), with periodic boundary conditions in all directions. The equilibration time was 
10 5 MD time steps. At this high temperature, only very weak chemical correlations develop 
among the particles. We use the standard velocity Verlet algorithm 95 ' 96,97 with a time step 
of 0.02, and apply the Nose-Hoover algorith m 95 ' 96 ' 97 for thermalization. 
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At time t = 0, the wall potentials are introduced, and the temperature is quenched to 
T = 1.1. This is done by rescaling the velocities, and by setting the temperature of the Nose- 
Hoover thermostat to the new temperature. Of course, in a real experiment, the temperature 
of a fluid confined in a small slit pore would be controlled via the thermal energy of the solid 
walls forming the pore. Therefore, no instantaneous quench (on picosecond or nanosecond 
time-scales) is possible. However, the structure formation occurring in a binary fluid with a 
finite quench rate is a complication that we disregard here. 

It is also relevant to discuss our procedure of introducing the walls together with the 
quench at time t = 0. In this case, the fluid is translationally invariant for t < 0, but 
loses this invariance in the z-direction for t > 0. However, we have found that the typical 
oscillatory density profiles near the walls ("layering") already develop during the first few 
MD time steps after the quench - see Fig. El For D = 5, we recognize 7 well-developed layers 
and there is no region of constant density in such an ultra-thin film. However, for D = 10, 
the region from z ~ 4 to z ~ 6 has an almost constant density n(z) = nj±{z)-\-n-B{z) ~ n = 1. 
For the D = 20 case, this constant density region covers about half of the film thickness, 
extending from z ~ 5 to z ~ 15. Note that the layer distance in Fig. OJis slightly less than a, 
although a coincides with the position of the first peak of the radial distribution functions 
in the bulk^l. We introduce the walls together with the quench at time t — to make 
the initial state of the quench (random distribution of A and B particles everywhere in the 
system, also close to the wall) comparable to that of the Ginzburg-Landau model (where we 
quench from a state at "infinite temperature"). 

It has been emphasized*^ that MD simulations constitute a method to explore hydro- 
dynamic phenomena, and are competitive with coarse-grained methods. In principle, this 
is only true for a microcanonical MD in the NVE ensemble where the energy E is strictly 
conserved. However, here we use the Nose- Hoover thermosta t 95 i 96 i 97 , i.e., we integrate the 
following equations of motion for the coordinates of the particles fi(t) (with r; = dfi/dt = Vi) : 



Here, Q is the fictitious mass of the thermostat, which was set to Q — 100. In the limit 



cm 



(18) 




(19) 
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Q = oo, we have ((t) = and then we recover the strict conservation of energy and 
momentum, on which the equations of hydrodynamics are based. For finite Q, we have 
{((t)) = 0. However, the fluctuating damping term disturbs hydrodynamics slightly. For 
this reason, in our earlier study of isothermal transport coefficients in the bul k 37 ' 38 , we have 
used strictly microcanonical runs. However, the ensemble of initial states was generated 
by Monte Carlo simulations in the semi-grand-canonical ensemble, ensuring thus a strict 
validity of all conservation laws in conjunction with averaging in the NVT ensemble. In 
the context of a thin fluid film confined in a slit between solid walls formed from vibrating 
atoms, neither momentum nor energy (of the fluid film) are conserved, and the walls do 
act as a thermostat. Of course, in a more realistic model, this thermostatting action of 
the walls applies only to those fluid particles which are close to one of the walls, and not 
on particles near the center of the film (which are only thermostatted indirectly via heat 
conduction). In situations far from local equilibrium (such as strongly sheared fluid a 97 ' 98 ), 
there is indeed a noticeable difference between the effects of wall thermostats and the Nose- 
Hoover thermostat. However, we do not expect any such problem here, since the time-scales 
for domain coarsening are much larger than the time-scales associated with heat conduction. 

In simulations of domain growth, one encounters the problem of large statistical fluc- 
tuations, and quantities such as the equal-time correlation function C(f,t) exhibit lack of 
self-averaging 99 . [Unlike the equilibrium case, C(f,t) explicitly depends on the time t after 
the quench.] Such quantities can only be sampled if a number of independent runs are per- 
formed and averaged over. All statistical quantities presented here are obtained as averages 
over three independent runs. 

IV. NUMERICAL RESULTS 

Let us begin with a discussion of the laterally averaged order parameter profiles, ipa, v (z, t) 
vs. z (see Fig. 0]). These are obtained from our MD simulations by averaging individual 
profiles for ip(x,y, z,t) vs. z along the x and y directions, and further averaging over inde- 
pendent runs. The morphology of these SDSD profiles consists of an A-rich wetting layer 
at the surface, followed by a depletion layer in A, etc. One can see that the order pa- 
rameter at the surface has already increased to a rather large value at early times, due to 
the preferential attraction of the A-particles to the walls. These A-particles are removed 



15 



from the adjacent regions in the interior of the film, resulting in local minima in if) av (z,t) 
(A-depletion layers). As time proceeds, these minima move towards the center, and also 
become more pronounced. In the D = 5 case, the SDSD waves coalesce rapidly and only a 
single minimum in the center is left by t = 640. In the D = 10 case, distinct SDSD waves 
are visible till t ~ 1000. At t — 2000 (see Fig. QJI), the waves have merged to give a layered 
structure with a single minimum. In both cases, the layered structure is transient and breaks 
up into a columnar structure which coarsens laterally (see Fig. Of course, even in this 
asymptotic state, the walls remain A-rich and the film center is B-rich. This behavior is rem- 
iniscent of the SDSD profiles seen in GL studies of this proble m 28 ' 31 ' 40 ' , and corresponding 
experiments 3 ^. In the present MD simulations, only a single depletion minimum is observed 
near each wall - statistical fluctuations of the local position of the boundaries between the 
depletion layers and adjacent enrichment layers wipe out any further systematic variation 
of the concentration profiles. 

It is also interesting to examine the evolution of the local order parameter ip av (0,t) at 
the surface z = (see Fig. Et) or z = D (which is analogous to Fig. Et)- Note that we have 
averaged the MD data over a layer of thickness Az = 1 to estimate ipa, v (0, t) (whereas for the 
calculation of ijj m {z,t) in Fig. ^\Az = 0.25 was used). The quantity ip av (0,t) rises rapidly, 
and reaches a maximum at about two decades, before it starts decreasing. The rapid rise 
is expected from the phenomenological theory (see Sec. II C). Of course, due to the lack 
of conservation of the local order parameter adjacent to the walls, there is an immediate 
response to the surface potential at the wall. For both D = 5 and D = 10, even runs 
up to t = 2 ■ 10 4 do not suffice to estimate the final values of ip av (0,t) clearly. This non- 
monotonic relaxation is a consequence of the structural rearrangement of the concentration 
inhomogeneities in the films. Figure shows corresponding data for ip av (0, t) vs. t from our 
recent simulations, using the GL model described in Sec. II. C 41 ". (Note that a logarithmic 
time-axis is chosen in Fig. Eb-) The behavior of the GL data is qualitatively similar to that 
of the MD results. However, a pronounced intermediate plateau is formed in the GL case, 
whereas the MD data only show a maximum - see the inset of Fig. Eb, which plots the data 
from Fig. on a logarithmic time-scale. The reason for this difference lies in the formation 
of a long-lived metastable layered state with pronounced A-rich layers in the GL case, which 
is not observed in the MD simulations. 

One way to further elucidate the morphological evolution in the MD simulations is to 
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look at snapshot pictures of the concentration in cross-section planes through the films (see 
Fig. |BJ). At t — 80, the distribution of the particles is rather random, but by t = 800 
the existence of domain structures is fairly evident. This interconnected structure coarsens 
(t = 1600), and ultimately breaks up into compact domains that connect the enrichment 
layers at both walls. (In the snapshot for D = 10, at t = 8000, only the enrichment layers 
are seen, but this observation is accidental. The slice shown cuts through a region free of 
columnar A-rich domains over the lateral scale L = 64, while other slices parallel to the one 
shown do cut through such domains. But we include this example in order to emphasize 
that strong fluctuations occur, not only from one run to the next run, but also in the course 
of the time evolution of individual runs. Thus, individual snapshot pictures give qualitative 
insight only.) 

Thus, Fig. El already provides evidence of the simultaneous presence of surface enrichment 
layers and lateral phase separation. A clear picture of lateral phase separation is obtained 
if we examine the concentration distribution in slices (of width a) centered at the L x 
L mid-plane at z = D/2 (see Fig. Ej). These pictures resemble snapshot pictures of 2-d 
spinodal decomposition, though for an off-critical composition. Of course, the concentration 
is conserved in a strictly 2-d system, whereas the concentration is not conserved in the shown 
L x L x a slice. As a matter of fact, it decreases systematically with increasing time, due to 
the progressive formation of A-rich surface enrichment layers. This is particularly evident 
in the late-time snapshots (t = 8000) for D = 10 and D = 20, respectively. 

In order to quantitatively characterize the lateral phase separation, we introduce the 
layer-wise correlation function: 



We also define a layer-wise length scale £(z, t) from the decay of this function with lateral 
distance p: 



In Fig. IHJ we plot the scaled layer-wise correlation function, C(p, z, t)/C(0, z, t) vs. p/£(z,t) 
at t — 8000, for D = 5, 10, 20 and different values of z. The surface (z = 0) is strongly 
enriched in the preferred component A - see Fig. 0] and Fig. |3| The corresponding corre- 
lation function measures small fluctuations about a strongly off-critical background. The 
correlation functions for the inner regions of the film do not scale either. (If there were 



C(p, z, t) = (V(0, z, t)1>(p, z,t))- (^(0, z, t)){ij(p, z, £)>. 



(20) 




(21) 
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scaling, all the data sets in Figs. |Hk-c would superimpose, as they do when similar plots are 
made in studies of bulk spinodal decomposition.) This is because the correlation function is 
a function of the off-criticalit y 100 ' , and different values of z are characterized by different 
average compositions - see the depth profiles for t = 8000 in Fig. (for D — 5) and Fig. EJi 
(for D = 10). 

In Fig. El we show the scaled correlation function in the film center for D = 5, 10, 20 
at different times. Again, there is no scaling of the data sets. This lack of scaling is 
expected, however, since the average volume fraction of A in the central region changes 
with time - see Fig. |3J For the case of D = 5, there is a reasonable superposition of the 
curves for t = 4000 and t = 8000. This is consistent with the observation that the average 
concentration at the center remains approximately unchanged over this time-regime - see 
Fig-lib- in the asymptotic regime, the system evolves via the lateral coarsening of columnar 
domains. Therefore, we expect the depth profiles ip aY {z,t) vs. z (as in Fig. 0j) to become 
independent of time at sufficiently large times. In this regime, we will recover dynamical 
scaling for the layer-wise correlation functions. 

In Fig. we plot the layer-wise length scale [£(z, t) vs. t] for D = 5, 10, 20 on a log-log 
plot, in order to check for possible power laws. At early times, no well-defined power law can 
be identified at all. This is not surprising as one does not expect a universal growth law to 
apply when the length scale is of the same order as the inter-particle distance. The gradual 
increase of the slope of d[ln£(t)]/d(lnt) is consistent with data from experimental and sim- 
ulational studies of spinodal decomposition in the bul k 21 ' 22 ' 2 ^ 24 ' 2 ^ 4 ^ 4 ^ 4 ?^ 1 ?^ 1 . Surprisingly, 
at later times, the MD data appear to be compatible with a power law with an effective 
exponent ~ 2/3. We do not see evidence for any of the other growth laws discussed in the 
context of fluids (see Sec. IIB) over an extended period of time. One might have expected 
that the LS evaporation-condensation mechanism or the droplet diffusion-coagulation mech- 
anism would dominate over some time-range, but this is not the case. As regards the Siggia 
tube-coarsening mechanism, the interconnected domain structures break up so early that 
hydrodynamic mechanisms can hardly become operative. 

It would be premature to claim that the log-log plots in Fig. El are evidence that the 
inertial mechanism [Eq. (J3j)] has been seen. According to theory, this mechanism should be 
visible only if the length scale £{t) ^> £- m = rf / (ti'Ya.b) • Fortunately, the material parameters 
which determine l- m are known for our model, as emphasized in Sec. III. Putting the numbers 
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in, we estimate that £- m ~ O(10 2 )! Such large values of £- m are compatible with studies of 
the late stages of domain growth in d = 3 using the lattice Boltzmann method 52 . 

One might then conclude that the effective power law £(t) oc t 2 ^ 3 seen in Fig. [HI] is only a 
transient phenomenon, and the power laws that are expected in this case (see Sec. IIB) will 
come into play at later times. An alternative possibility is that novel growth laws arise due 
to the interplay of wetting kinetics and lateral phase separation in the "bulk" of the film. 
However, we note that our results have a striking qualitative similarity to the Brownian- 
dynamics results of Farrell and Valls^i. These authors studied phase separation in strictly 
2-d fluid mixtures, and found a rapid crossover to a power law £{t) oc v' with ~ 2/3. 

Obviously, more work with both simulations and theory is needed to resolve the nature 
of applicable growth laws. However, this cannot be done by simply running our simulations 
longer. This is because the condition £(z, t) <C L is needed to ensure that finite-size effects 
in the lateral direction are negligible. Further, the condition £(z, < L is also needed 
to provide a reasonable self-averaging of C(p,z,t) [defined in Eq. (J20|) ]. One can divide 
the system laterally into independent blocks of linear size £(z,t) to judge the error in the 
estimation of C(p,z,t). Therefore, the relative error is less for D = 5 (where L = 128) 
than for D = 10 and 20 (where L = 64), and it increases when £(z,t) increases. Thus, 
the irregularities in £(z,t) for D = 10,20 when t > 1000 are probably due to insufficient 
statistics (as only three independent runs were made). 

V. SUMMARY AND DISCUSSION 

Let us conclude this paper with a summary and discussion of our results. Here, we have 
presented comprehensive results from molecular dynamics (MD) simulations of surface- 
directed spinodal decomposition (SDSD). We have used a simple model system, namely a 
symmetric binary Lennard- Jones (LJ) mixture, confined between identical flat and struc- 
tureless parallel walls which preferentially attract the A-particles. Only very thin films are 
accessible - the distance between the origins of the wall potentials was (D + l) = 6, 11, 21 in 
units of the LJ parameter. Further, the finite size of the lateral linear dimension L (L = 128 
for D = 5, and L = 64 for D = 10, 20) constrains our work to the early and intermediate 
stages of domain growth. In this regime, the characteristic length scale of lateral phase 
separation £(z, t) has grown by approximately one decade. Note that we have also restricted 
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attention to deep quenches, much below the critical temperature of phase separation, but 
above the triple-point temperature, so crystallization is not an issue in our study. At the 
chosen temperature (T = 1.1, i.e., T/T c ~ 0.67), the bulk phase separation occurs between 
almost pure A and B fluids. The interfaces are locally sharp (with a correlation length 
£ ~ 1), and the time-scale for structural relaxation in the fluids is manageable for MD work. 
The shear viscosity has been estimated previously^ to be 77 ~ 7 at T = 1.1, in the standard 
LJ units. Thus, the advantages of the present approach are as follows: (a) all material 
parameters of the model are explicitly known; (b) at small scales, a qualitatively reasonable 
description of fluid structure is ensured; and (c) long-range hydrodynamic interactions (re- 
sulting from the conservation laws in fluid dynamics) are automatically included, although 
only a short-range LJ potential is chosen to model the interaction among the point particles. 

We use this model to elucidate all the main characteristics of SDSD. The surfaces become 
the origin of SDSD waves, which consist of alternating enrichment and depletion layers of 
the preferred component A. These waves coalesce in the central region of the film, giving rise 
to a layered structure - see the t = 800 profile in Fig. 0]d and the t = 4000 profile in Fig. EJi. 
This layered state subsequently breaks up into a columnar structure that coarsens laterally 
- see the cross-sections in Figs. |B] and Therefore, the local concentration of A-particles at 
the walls grows rapidly at first, resulting in rather large values at early times, followed by a 
decrease at later times (Fig. EJ). 

In the initial stages of phase separation, domain growth in the film interior resembles 
spinodal decomposition in bulk mixtures, where a bicontinuous percolating structure forms 
at compositions near the critical concentration. During this stage, £(z, t) grows only rather 
slowly (Fig. irUj) . In this regime, the concentration of A in the film center decreases, due 
to the growth in thickness of the surface enrichment layers. Therefore, the percolating 
structure breaks up into separated A-rich droplets which have a cylindrical shape, with 
height of order D in the z-direction and radius of order £(z,t). However, these droplets are 
connected through the A-rich enrichment layers at the walls. The observation that, in this 
droplet growth stage, the local concentration at the surface decreases can be understood as 
follows. For the chosen parameters, the B-rich phase exhibits incomplete wetting of A at the 
walls - for complete wetting, no overshoot of ij) av {z = 0, t) vs. t would be expected in Fig. |SJ 
A remarkable feature of our results is that domain growth is compatible with the inertial 
growth law, £(t) oc t 2 ^ 3 , during the droplet growth stage (Fig. fTTT^ . This is reminiscent 
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of studies^ of phase separation in strictly 2-d fluids, where Langevin equations including 
hydrodynamic interactions were simulated. 

It is clear that an extension of the brute force MD approach to much larger linear di- 
mensions (needed at later times) requires prohibitively large amounts of computer time. On 
the other hand, our MD study does reach mesoscopic length scales, which are significantly 
larger than inter-particle distances. This suggests that our MD study should be supple- 
mented by Langevin studies of the Ginzburg-Landau (GL) models described in Sec. II. C, so 
as to enable a simulation extending from microscopic to macroscopic scales. Let us briefly 
discuss a comparison for the early stages of SDSD in thin films, where the hydrodynamic 
interactions can be disregarded. Figure ^2 is analogous to Fig. but is taken from a GL 
simulation of model B with appropriate boundary conditions 41 . The details of this simula- 
tion are discussed at the end of Sec. II. C. The qualitative similarity of the evolution of the 
depth profiles in Figs. 0] and ^2 is striking. For the GL simulations, the spatial degrees of 
freedom were rather coarsely discretized (with Ax = 1), and hence the small-scale structure 
close to the walls cannot be resolved. Apart from this difference, the GL profiles are in good 
agreement with those obtained from the MD simulation for short times, if we equate the 
time units as £ql — 160 £md- Recall that the natural time-scale in fluids is given in terms 
of the structural relaxation time^., and the latter is of the same order as the shear viscosity, 
which is r\ ~ 7 LJ units in this case. Since the MD time unit £md = — l/v^48 ~ 1/7 L J 
units, we conclude that the natural fluid time-scale is about 50 MD time units. 

The qualitative behavior of the GL results 41 - (formation of surface enrichment layers 
and a metastable layered state — ► break up due to lateral phase separation — > coarsening 
of columnar structures) is similar to the results of the present MD study. However, in 
quantitative respects, the intermediate stages of coarsening are rather different for the 
GL model. When the layered state breaks up into a laterally inhomogeneous state in the 
GL model, one often encounters a period of time where i stays roughly constant, or even 
decreases. We do not observe such a transient behavior in the MD runs - as a matter of fact, 
the layered state does not survive for any appreciable period of time. Further, asymptotic 
domain growth in the GL model is compatible with an i(t) oc t 1 / 3 law, as predicted by 
the Lifshitz-Slyozov theory. On the other hand, our MD results are consistent with a 
growth law £(t) oc t 2 ^ 3 . Clearly, a GL study of SDSD in thin films, which incorporates 
hydrodynamic interactions, would be very desirable. Further, a detailed comparison of the 
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present MD results with corresponding lattice Boltzmann studies should be worthwhile, 
but is left to future studies. Finally, we hope that the present work will further stimulate 
experimental studies of phase separation in thin films. 
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D=20 



t=80 




FIG. 1: Snapshot pictures of surface-directed spinodal decomposition (SDSD) in a binary (AB) 
Lennard-Jones (LJ) mixture, which is confined in an L x L x D thin-film geometry, with L = 64, 
D = 20. (All lengths are measured in units of the LJ diameter.) Periodic boundary conditions 
were applied in the lateral directions, while the impenetrable Lx L surfaces (representing the walls 
of a slit pore) attract the A-particles. The initial condition for this run consisted of a random 
mixture of equal amounts of A and B, corresponding to a critical quench. Time is also measured 
in dimensionless LJ units - see Sec. Ill, where further details of the simulation are specified. The 
A-particles are marked black, and the B-particles are marked gray. The system quickly develops 
concentration inhomogeneities (t = 800), in particular A-rich layers form rapidly at the walls 
(t = 800, t = 1600). The late stages of phase separation are characterized by the lateral coarsening 
of columnar structures (t = 8000). 
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FIG. 2: Schematic phase diagrams (a,b) and corresponding states (c) of a symmetric AB mixture 
in a thin film of thickness D. The film is symmetric, viz., both walls attract the A-particles with 
the same strength. We emphasize that the wetting transition only occurs in the limit D — > oo. For 
finite D, the transition of the walls from nonwet (or partially wet) to wet (or completely wet) is 
rounded into a smooth gradual change. This transition is of second order in (a), while (b) refers to 
a first-order wetting transition. In (b), a prewetting transition line exists in the one-phase region, 
with one end being a prewetting critical point at high temperatures. The other end of this line is at 
the wetting transition temperature T w , at the coexistence curve that separates the two-phase region 
from the one-phase region. Note that the critical concentration of a symmetric binary mixture is 
x™* = 0.5 in the bulk, but is shifted to a larger value xa in the thin film. Further, the critical 
temperature of the film typically is lower than in the bulk, T C {D) < T c (oo) = T c . For the case 
of first-order wetting and large enough D, a thin-film analog of the prewetting transition exists, 
as evidenced by the thin-film critical point at the left side of the phase diagram. When the thin 
enrichment layer segregation meets the lateral segregation of the "thick" film, a thin-film triple 
point occurs at a temperature close to T w . For thin films, this triple point and the left critical 
point may merge and annihilate each other, and then the corresponding phase diagram is similar 
to that in (a). In (c), we provide schematic pictures of the thin-film states in the case of lateral 
phase separation. 
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FIG. 3: Profiles of the total density n(z) = n^{z) + ns(z) vs. z for (a) D = 5, (b) D = 10, and 
(c) D = 20. We show data for four different times t (in units of to), as indicated. The vertical 
line in each frame indicates the mid-plane of the film at z = D/2. The wall potentials diverge at 
z = —1/2 and z = D + 1/2 [see Eq. (|14|)]. so the particles can range over a distance (D + 1) in the 
z-direction. All lengths are measured in units of a, and hence are dimensionless. 
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FIG. 4: Laterally averaged order parameter profiles, ip ay (z,t) vs. z, for films of thickness (a) 
D = 5 at early times (t = 80, 160, 320, 640); (b) L> = 5 at late times (t = 800, 2000, 4000, 8000); (c) 
D = 10 at early times (t = 80, 160,320,640); (d) D = 10 at late times (t = 800,2000,4000,8000). 
The symbol usage is the same for (a),(c) as well as for (b),(d). 
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FIG. 5: (a) Time-dependence of the local order parameter ^avCO, t) at the surface z = 0. We show 
data for D = 5 and D = 10. Note that we have averaged the MD data over a layer of thickness 
Az = 1 to estimate ^ av (0,t). (b) Time-dependence of Vw(0, t), obtained from the Ginzburg- 
Landau (GL) simulations described in Sec. II. G— . The GL data was obtained as an average over 
5 independent runs with L = 256. We plot the data on a log-linear scale. The inset shows the MD 
data from (a), also on a log- linear scale. 
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FIG. 6: Snapshots of the concentration distribution in cross-section slices (of linear dimensions 
a x L x D) through the films, for (a) D = 5, (b) D = 10, and (c) D = 20. The cross-section 
was centered at x = L/2. The A-particles are marked black, while the B-particles are not shown. 
These pictures correspond to the times t = 80, 800, 1600, 8000 in all cases. 
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FIG. 7: Snapshots of the concentration distribution in cross-section slices (of linear dimensions 
Lx Lxa) centered at the plane z = D/2 of the films. We show pictures for (a) D = 5, (b) D = 10, 
and (c) D = 20. The A-particles are marked black, while the B-particles are not shown. Note that 
L = 128 for D = 5, but L = 64 for D = 10, 20. 
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FIG. 8: Plot of normalized correlation function C(p, z,t)/C(0, z,t) vs. p/£(z,t) at t = 8000 for 
(a) D = 5, (b) .D = 10, and (c) I? = 20. We show data for several values of z (distance from the 
left wall), as indicated in the figure. The time t = 8000 is chosen such that it corresponds to the 
later stages of coarsening. 
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FIG. 9: Plot of normalized correlation function C(p, z,t)/C(0, z,t) vs. p/£(z,t) for z = D/2, and 
(a) D = 5, (b) D = 10, and (c) D = 20. We show data for times t = 80, 800, 4000, 8000 in all cases, 
as indicated by the different symbols. 
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FIG. 10: Time-dependence of layer- wise length scale £(z,t), plotted on a log- log scale. We show- 
data for different values of z, and (a) D = 5, (b) D = 10, (c) D = 20. The straight lines have a 
slope of 2/3. 



36 



(a) D=5 (GL) (b) D=5 (GL) 




012345 012345 

z z 

(c) D=10(GL) (d) D=10(GL) 




02468 10 02468 10 

z z 

FIG. 11: Laterally averaged order parameter profiles, ip ecv {z,t) vs. z, obtained from the GL 
simulations described in Sec. II.G^k The GL data was obtained as an average over 5 independent 
runs with L = 256. We show data for films of thickness (a) D = 5 at early times (t = 0.5, 1, 2, 4); 
(b) D = 5 at late times (t = 10,100,1000,20000); (c) D = 10 at early times (t = 0.5,1,2,4); (d) 
D = 10 at late times (t = 10, 100, 1000, 20000). The symbol usage is the same for (a),(c) as well as 
for (b),(d). 
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